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We present two improvements to our previous dynamical overlap HMC algorithm. We introduce 
a new method of differentiating the eigenvectors of the Kernel operator, which removes an in- 
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without pseudo-fermions, we are able to increase the rate of topological tunnelUng by a factor of 
more than ten, reducing the auto-correlation. 
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1. Introduction 

The overlap operator ^ is the only known lattice Dirac operator with an exact lattice chiral 
symmetry Since chiral symmetry is important for many low energy observables, it is desirable 
to use overlap fermions. However, using a full QCD simulation with overlap fermions presents 
a number of algorithmic challenges. In this paper, we address, and present solutions for, two 
outstanding algorithmic issues. 

The problem of the Dirac 5-function in the fermionic force when changing topological charge 
has been resolved by using a transmission/reflection algorithm, similar to the case of a classical 
mechanics particle approaching a potential wall. The original formulation [^, ^ has subsequently 
been improved in to maximise the rate of topological charge change for a given action jump. 
Our method is described in these references. 

To differentiate the overlap operator in the Hybrid Monte Carlo (HMC) molecular dynamics 
(MD), it is necessary to differentiate the eigenvectors and eigenvalues of a sparse matrix. Previous 
methods have led to instabilities when there are degenerate eigenvalues. In section || we discuss a 
new method which avoids these instabilities 

The topological auto-correlation depends on the rate of topological activity. To have high rate 
of topological index change, it is necessary to reduce the discontinuity in the action at this point. 
In section ^, we outline a new method for reducing this action jump [^. 

Our overlap operator is defined as 

D=(l + At) + 75(l-iU)sign(!2). (1-1) 
In our tests, Q is the standard Wilson operator, 

e = 75 [5xy - ((1 - 7^)^/^: W5,,.v+M + (1 + Y^W^X - pi)5y,,-y,)] , (1.2) 

with K = 0.2 and, in section |3[ two levels of stout smearing [||] at parameter 0.1. Our numerical 
tests are performed on 8^16 lattices at a lattice spacing of about O.I5fm (measured using ro), with 
quark masses /i = 0.03,0.04,0.05, corresponding to pion masses in the range 500 — lOOOMeV. 

2. Eigenvector mixing 

2.1 Differentiating Eigenvectors 

The eigenvalues A,- and eigenvectors ) of a matrix Q are defined by 

QWi) = UWi)- (2.1) 
After a small change to the matrix, 5Q, the new eigenvalue equation is 

{Q + 8QM)=XM)- (2.2) 
We can expand the new eigenvectors in terms of the old basis 

\w'i) = Wi)+ L(cos % - 1 ) I V^,-) + e'^'' sin dtj (2.3) 
j 
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We assume that only one of the mixing angles dij is large, so it is unimportant that ly/J) is not 
normalised^ . The mixing angles dij and ^ij are 



tan 2 a- 



l^{^|fi\5QWj){^i^j\5Q\^ri) 
''Xj-Xi + {WjmWj)-{WimWi) 
l{Vj\5Q\Vi) 



V {¥i\SQ\Vj)' 
Expanding 6 and in r/{Xi — Ay) gives to lowest order 

[Wj\5Q\¥i) _ 1 



(2.4) 



Xj — Xi 



\\ifi){Wi\)^QWi), 



(2.5) 



which agrees with other methods (see, for example, [g]). Equation ( |2.5| ) breaks down when A,- — Xj 
is small. In this situation, it is necessary to use the exact expressions for the mixing angles (2.4). 
Defining (x) as the eight SU(3) generators on one link, we can write 



'■ij > 



5Qij = {Yi\5Q\wj) = r7i';,{x)al 
where n represents the MD momentum, and the vectors are given by 



(2.6) 



a. 



n,x,ll 



M¥i\r5[{l-r^iK{x)u^{x)5y,,,+^-{l + r^)ul{x)T|,{x^ iu) 



The generalisation to the smeared operator is trivial The NAC force, F^'^^ , constructed for the 
overlap operator in |^, is proportional to a,y, and a function of only na and the gauge field, U . 
We can construct a reversible algorithm by combining forward and backward half-steps (using a 
quick and reversible iterative procedure for the backward step). The force is not area conserving, 
but it is possible to use a non-area conserving force by calculating the Jacobian, J, and including 
log/ in the HMC accept/reject step (see |^ for an example). J can be written as 



3F 



n+Aiikiaijau 



It is now a simple task to rewrite a,y in terms of an orthonormal complete basis a' , and use 



det 



:det[l+A'] 



(2.8) 



(2.9) 



to calculate J. log / scales with T^, and does not affect the HMC acceptance rate. 
2.2 Numerical results 

In figure |l], we compare the fermionic forces for the old algorithm and the NAC algorithm 
across a typical HMC trajectory (using overlap fermions with no smearing). It is clear that the old 
method gives an unstable force, which cannot be used, while the NAC force is stable. In all our test 
trajectories, we have not observed a log / larger than 0.3. 



A more general expression can easily be constructed should this assumption break down. 
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Figure 1: Comparison of the trace of the square of the fermionic forces for the proposed and old algorithms 
with T = 0.016 on one of the pL = 0.05 trajectories. 

3. Determinant factorisation 



3.1 Introduction 

The topological auto-correlation is related to the rate of topological activity. We need to be 
able to measure this rate, from which we can estimate the auto-correlation. Since instanton anti- 
instanton pairs are difficult to observe, our best way of measuring the auto-correlation is from 
topological charge changes. The probability of transmission scales as mm{\,e^^) [||], where AS', 
the action discontinuity at the topological sector boundary, scales as /i^^. Therefore it is necessary 
to reduce AS. 

However, the determinant of the actual Dirac operator, logdet(DD^), scales as log }JL rather than 



-2 



p. The difference is caused the different functional forms of the pseudo-fermion estimate and 



the actual determinant. It has already been shown that adding additional pseudo-fermions reduces 
AS considerably [1^, 11 1, and this (with one additional pseudo-fermion) is our algorithm A, which 
we test against our new methods. We propose [Q] factorising the determinant into a large continuous 
part, which can be treated with pseudo-fermions, and a discontinuous determinant small enough to 
be calculated exactly. 

3.2 Algorithm C 

Algorithm C factorises the determinant using 



det 



75 



1+M 
1-M 



+e(e) 



=det 



1-At 

= det[Di]det[D2] 



det 



Iv^;)(£(^y)-£(%) 



(3.1) 
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£ is a continuous approximation of the sign function, only differing for the smallest n eigenvalues 
below a cutoff A. In our tests, we used a Zolotarev rational approximation. We simulate det[Di] us- 
ing pseudo-fermions; det[D2] is calculated by standard methods. We add — logdet[D2] to the HMC 
action. To maintain a high HMC acceptance we calculate logdet[D2]'s force (unlike the proposal 



in 1 12]). Differentiating logdet[D2] is straight-forward using the methods of section Algorithm 



C requires an inversion of Di for each eigenvalue projected. While deflation methods [ |13| ] reduce 
this cost, it could lead to difficulties on larger volumes.^ 

3.3 Algorithm F 

To avoid this additional cost, algorithm F projects out just one vector, \a), which is equal to 
the smallest eigenvector at the moment of crossing. Thus, we write 



det[l + 75e(e)] =det[l + 75e(e)(l-|a) (a|)]det 

= det[D3]det[D4] (3.2) 



1 H 1 > , MX Ys^iQ) \a) (a| 



To ensure D3 is both continuous and practical, \a) must satisfy la) = |i//,) at A,- = and = 

for Xf > A^ (where A^ is some suitable eigenvalue cut-off). For the algorithm to remain stable 
at larger lattice volumes \a) and its differential must be continuous, d/dXj \a) must be sufficiently 
small, the eigenvalues and eigenvectors must be differentiated using the results of section ^ and the 
approximate overlap operator D3 must not have any exceptional configurations. We construct \a) 
from n eigenvectors using 

k)=A,lwj^ + /3i|v.i)j^ + -, (3.3) 
where F is a constant vector used to fix the relative phase of the eigenvectors, and currently we use 

tan^A- = A,^A^("-i)/(a2 - Af)2n^y ' (3-4) 

although we are still searching for the best function to use for j8 . Algorithm C has shown stable 
MD with high HMC acceptance and good reversibility on our test 8^16 ensembles.^ Algorithm F 
has encountered large differentials of j3 with respect to the eigenvalues, which may cause problems 
on larger volumes, although the HMC acceptance rate is still acceptable on our lattices. We do not 
expect any scaling of the action jump with the lattice volume. 

3.4 Numerical results 

Figure ^ and table |I| give preliminary results for the action jump on our 8^16 test configu- 
rations. A5 is much smaller for algorithms C and F than for algorithm A, and the distribution is 



^We cannot decrease A without increasing tiie force; tiierefore tiie number of projected eigenvalues will increase as 
the volume increases. 

^ A variant of algorithm C, not using our method of differentiating the eigenvectors and not treating the small matrix 
in the molecular dynamics, has encountered problems on 12-' 24 volumes ||l4]]. We suspect that this was because the 
range of the approximate sign function was set to small, and a less efficient method was used when differentiating 
eigenvectors. 
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Figure 2: The distribution of the action discontinuity A5 for algorithms C and F (left), and A, C and F 
(right), at masses pi = 0.03,0.04 and 0.05. 
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Table 1: The mean values of AS for algorithms A,C and F at three quark masses pi (left) and the standard 
deviations (right). 



much narrower. The results for both algorithms C and F give a high transmission rate for all /i. 
We see little variation with the quark mass with these methods, unlike algorithm A. This suggests 
transmission may be possible with much smaller /i. 

4. Conclusion 

We have developed two additions to the overlap HMC algorithm designed to significantly 
improve simulations at large volume and small fermion mass. By employing a non area conserving 
algorithm to differentiate the eigenvectors and eigenvalues of the Dirac operator, we remove various 
instabilities that otherwise are encountered in the fermionic force. By factorising the fermion 
determinant, we decrease the action jump at the topological sector boundary considerably (we have 
observed transmission rates of up to 80%), reducing the auto-correlation by an order of magnitude. 
Thus, with dynamical overlap fermions, it is possible to accurately sample all topological sectors 
efficiently. It remains an open question whether other lattice formulations will, at small lattice 
spacing, suffer from the large auto-correlations which we have now avoided. 
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